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I. INTRODUCTION 


Multiphase flows are often present in numerous 


Abstract — This work presents a methodology for reconstructing multiphase 
flow electrical capacitive tomography (ECT) images, using a particle swarm 
optimization (PSO) algorithm in the parallel processing paradigm. Intended 
is to improve the efficiency of the inverse problem algorithm in ECT, 
increasing the resolution of the reconstructed images, without necessarily 
increasing the processing time of these reconstruction technique. A 
limitation found is that, for inverse problem-type reconstruction techniques 
for ECT, the response of the sensor system is non-linear and, therefore, the 
processing time grows faster than any increase in resolution, imposing a 
high computational cost. For real-time applications, the first contribution is 
the removal of unnecessary processing from the usual code; the second is 
the creation of a new PSO algorithm for image reconstruction that is more 
efficient than normal. The new parallel processing routine present the 
physical principles of ECT, the heuristic algorithms used in the 
reconstruction process and the main concepts for parallel computing. 


some advantages, such as i) been non-radioactive, ii) non- 
invasive and non-intrusive, iii) cost effective, iv) non 
susceptible to adverse temperature and pressure conditions 


industrial processes, particularly in food, chemical and 
petroleum industries and also in energy plants, among 
others (Thorn, Johansen and Hjertaker, 2013; Zainal- 
Mokhtar and Mohamad-Saleh, 2013; Cui et al, 2014; Mei et 
al, 2016). In the majority of these studied flows the 
dielectric nature of matter dominates, and consequently, the 
use of electrical capacitance tomography (ECT) is 
reinforced as one of the most advantageous tomographic 


and v) relative fast response. 


This technique is based on measurement of capacitance 
changes in dielectric distribution of the material present in 
the multiphase flow. Electrodes are placed around the flow 
pipe, usually 8 to 16, and the capacitance values across each 
pair are measured. Those capacitance measurements are 
then applied to an appropriate reconstruction algorithm 


method to obtain images from such industrial processes. In which produces an image of the dielectric spatial 
comparison to other tomographic techniques, ECT offers 
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distribution and therefore the material across the sensor 
(Yang and Peng, 2003). 


The ECT sensor is composed by a number of electrodes 
that can be made of a cooper sheet placed around the flow 
pipe in a non-parallel way. Besides the capacitive sensor, 
there are also i) an electronic system of transduction to read 
a voltage signal which is proportional to the capacitance and 
between each pair of electrodes, and ii) a control and 
acquisition system that coordinates the multi-electrode 
measurements and converts the analog signal into the data 
to be used by the image reconstruction algorithm. 


Il. METHOD 
2.1 Forward and inverse problem of ECT 


Due to its specific features, ECT is considered as a soft- 
field tomography. Therefore, the material to be imaged 
modifies the electric sensing field. Such effect is highly 
non-linear and the image reconstruction process becomes 
more complex in comparison to conventional x-ray 
tomography (Belo, 1999). 


Image reconstruction using ECT is a two steps task: 
firstly, the forward problem should be solved, in other 
words, it is the determination, using the excitation data, of 
the electric potentials inside the domain and the respective 
response along the contour. Secondly, the inverse problem 
should then be solved, which is the determination of the 
dielectric distribution into the domain from the relation 
between excitation data and the boundary response. 


Forward Problem Solution 


The solution of the forward problem associated with 
ECT can be defined as the determination of the electric 
potential, and thus the capacitance between the electrodes, 
associated to a given dielectric distribution over the area of 
across section of the pipe in study. 


If there are no free charges within the section of the 
imaged region, the relationship between the capacitance and 
permittivity distribution is modeled by Equation (1): 


Vle,e(x, yVo(x, y)| =0 (1) 


where (x, y) E€ N represents the space domain, £ọ is the 
vacuum absolute permittivity, e(x,y) is the relative 
permittivity distribution of the material inside the imaged 
region, (x,y) is the electric potential and F is the closed 
curve enclosing the electrodes surrounding the region. 


The boundary conditions are defined when an electrode 
is excited with a potential g = Vo (source electrode) and all 
others are kept on ground level (g = 0). The measuring 
process is then multiplexed for each grounded electrode 
(sensing electrode) in order to obtain the respective 
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capacitance. Therefore, the potential distribution is 
expressed as 


a Vo, (x,y) E€ Ni 
Py) = ety EN, k #i (2) 
and for the capacitance we have 
Æ = 
C= Sr Eœ yW ol, ydan (3) 


As long as there is not a general analytic solution for 
all configurations of permittivity distribution and electrodes 
excitation conditions, only particular solutions can be 
achieved for specific configurations. For example, 
numerical routines using finite element method (FEM) have 
been developed to determine the capacitance values from 
the permittivity distribution and excitation profile of 
electrodes surrounding a pipe (Souza, 2009). 


FEM was chosen due to its versatility and robustness 
in dealing with complex geometries and inhomogeneous 
media and this is also the method implemented here to solve 
the forward problem. 


Inverse Problem Solution 


Solution of the inverse problem associated to ECT can 
be defined as the determination of the dielectric distribution 
over the area of a cross section of the pipe under analysis 
which is related to the capacitance values obtained by the 
sensor. 


One method of accomplishing this is to compare 
estimated values obtained from a numerical model with 
experimental measurements from a capacitance sensor. 
Thus, the traditional tomography problem is moved to the 
minimization of an error functional. Such functional should 
be able to reflect the discrepancies between the 
experimental measurements values and numerically 
calculated values corresponding to changes in permittivity 
distribution inside the region to be imaged. 


The optimization process consists in a search algorithm 
for the dielectric distribution that globally minimizes the 
error functional (Smolik, 2010). At initialization, the 
forward problem is solved for a given initial distribution 
which results in a numerical answer to the capacitance value 
that gives a new value to the error functional. If such value 
is smaller than the earlier value, then the former distribution 
is updated and this process will be iterated many times until 
a global minimum for the error functional is found. When 
the process ends the ultimate numerical result obtained 
should correspond to the distribution that most resembles to 
the actual permittivity distribution inside the sensor, in other 
words, an image of the sensor cross section is achieved. 


When choosing a reconstruction algorithm, the main 
considerations that we should take into account are 
computational effort, speed and accuracy. Iterative 
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algorithms can increase accuracy but, on the other hand, 
may slow the rebuilding process. Image reconstruction 
process using iterative approach is time consuming, because 
it has to estimate the capacitance value by solving the 
forward problem many times, in an iterated fashion, thus 
becoming a low-speed reconstruction technique, one that 
cannot be applied to systems that require real-time images. 


2.2 RECONSTRUCTION ALGORITHM 


Unlike direct algorithms, the iterative algorithms are 
formulated in terms of an optimization problem, 
characterized by iterated attempts to minimize an error 
function between the capacitance values obtained by the 
theoretical solution of the forward problem and 
experimental capacitance values obtained by the system 
sensor. It is assumed as the solution of the inverse problem, 
the dielectric configuration that can minimize the error 
functional (Li, 2015). 


Particle Swarm Optimization 


Particle Swarm Optimization (PSO) is a stochastic 
computational technique based on population dynamics. 
PSO has emerged from experiences developed by Kennedy 
and Eberhart (1995) with algorithms that model a social 
behavior from a set of individuals. 


Likewise, other collective intelligence approaches, PSO 
is based on the skills of a population of individuals which 
are able to interact with each other and also the 
environment. Relying on the capabilities of self-evaluation, 
comparison and imitation, individuals are able to deal with 
a number of possible situations presented by the 
environment and therefore, global behaviors emerge from 
these interactions. The algorithm developed by Kennedy 
and Ebehart (1995) seeks the optimization of a fitness 
function using information exchanges among individuals 
(particles) and the whole population (swarm). 


In order to achieve the optimal solution, each particle in 
PSO is treated as a point in R” space and it represents a 
potential solution of the forward problem. Particle position 
is adjusted according to its own experience as well as group 
experience. Equation (4) corresponds to the sum of three 
distinct terms: the first one refers to the particle inertia; the 
second is a cognitive term related to particle individual 
learning of the best position it had already reached; the third 
is a social term that represents the experience exchange 
among all particles. In each iteration, the particle position is 
updated according to Equation (5) which considers its 
current position as well as a displacement given by the 
velocity due to iterating process. 


Vi = WY; + C11 (Ppest — Xi) +C27%2(Goest — Xi) (4) 


Xi = Xj + Vi (5) 
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Where v; is the particle i current velocity, wis an inertial 
weight that balances global and local exploration, c; and c2 
are behavior coefficients, r; and r2 are random numbers 
between Oand J, Poesr is the best position already reached by 
the particle and goes is the best position found by the swarm. 


Velocity update of particles depends on parameters that 
should be adjusted for each problem to be optimized, 
namely the inertial weight, cognitive and social terms. 
Inertial weight w allows broadness on the exploration space: 
high values improve global exploration. To the contrary, 
small values favor local exploration. In this work, the 
inertial weight has an update procedure (Equation 6) 
identical to what is described in Eberhart and Shi (2001), 
where w is adjusted linearly in the interval 0.4 to 0.9. Shi 
and Eberhart (1998) suggested to keep c; = c2 = 2.0 in order 
to balance social and cognitive behavior of the particle. 


(k-1) 
m 


(6) 


Where k is the iteration number, m is the maximum 
iteration number, Wmax and Wmin are the maximum and 
minimum weight, respectively. The steps for implementing 
PSO algorithm are described in the Pseudocode 1 (Eberhart 
and Shi, 2001): 


W = Wmax — (Wmax > Wmin) 


Pseudocode 1. PSO algorithm 


1 initialize constants M, N, C1, C2, W; 
2 randomly initialize n positions and velocities in the space R”; 
3 fork — 1tomdo 
4 fori — 1 ton do 
5 solve forward problem for position i; 
6 evaluate fitness f(x;); 
7 end for 
8 fori — 1 ton do 

seri i 
9 iff < f , then 

i _ gi 
=i f= fk 
=e 

a Prest = Xx 
12 end if 

spri 
13 iff < f ‚then 

_ ¢i 
US f= fk 
_ vi 

1> Ibest = Xx 
16 end if 
17 end for 
18 update velocity according to Equation (4); 
19 update position according to Equation (5); 
20 update w according to Equation (6); 
21 end for 


Result: assume gpes:as the problem solution; 


Parallel algorithm implementation of PSO can be 
described as a master-slave paradigm (Schutte et al, 2004). 
The master processor (MP) creates a set of random initial 
positions in the space R”, partitions this set into subsets and 
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sends them to have its fitness function analyzed by the slave 
processors (SPs). All decision processes are carried out by 
the MP, for instance, position and velocity update as well as 
the algorithm convergence control. SPs solve the forward 
problem and compare their results to the capacitance values 
obtained experimentally. Thus, for each position 
configuration the fitness function is evaluated and the 
resulting value is returned to the MP. MP and SPs tasks are 
described in Pseudocode 2. 


Pseudocode 2. Parallel PSO algorithm 


1 if my_id = 0 then // I am the master 

2 initialize constants itermax, N, C1, C2, W; 

3 randomly initialize n positions and velocities in the space R”; 
a endif 

s fork — 1 to iterma do 

6 if my_id = 0 then = //I am the master 

7 for target — 1 tO Nprocessors do 

8 send position x; of m particles to each slave 

9 


end for 
10 for source — 1 tO Nprocessors dO 
11 receive position Xi of m particles from each slave; 
12 receive fitness f(x;) of m particles from each slave; 
13 end for 
14 perform barrier synchronization for all fitness evaluation results; 
15 fori — 1 to n do 
16 if f,<f,, then 
17 Cah 
a Poest = Xk 
19 end if 
20 if f< f then 
a mi 
as g best — X k 
23 end if 
24 end for 
25 update velocity according to Equation (4); 
26 update position according to Equation (5); 
27 update w according to Equation (6); 
28 else // I am a slave 
29 receive position x; of m particles from the master; 
30 fori — 1 to m do 
31 solve forward problem for each Xi; 
32 evaluate fitness f(x;); 
33 end for 
34 send fitness f(x;) of m particles back to the master; 
35 end if 
3 end for 


Result: Assume gpes: as the problem solution; 


The initial steps of serial and parallel PSO algorithms 
are the same. Once the initialization step has been carried 
out by the MP, a package of m particles is sent to the SPs 
where the fitness function for each particle is evaluated. 


Information exchange between MP and SPs is 
implemented using Message Passing Interface (MPI), a 
peer-to-peer communication library, as described in Singh 
(2012). 


Parallel algorithm application aims to reduce processing 
time span from capacitance measurements (from sensor) to 
image reconstruction in comparison to serial algorithm. 
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Such reduction is attained by the distribution of forward 
problem solving to many processors. 


2.3 EVALUATION OF ALGORITHMS BY 
SIMULATION 


Algorithm process 


To analyze the spatial resolution obtained by the PSO 
algorithm applied to the ECT, we numerically simulated the 
response of a capacitive sensor with eight electrodes, with 3 
different distributions of a two-phase mixture (oil-water): A 
core flow distribution, a homogeneous distribution of oil 
water containing bubbles and finally a stratified 
distribution. 


Capacitance values were obtained from the direct 
problem solution for simulated values of permittivity, where 
they were contaminated with numerical noise in order to 
observe the impact of experimental error in the quality of 
the reconstructed image. The relative permittivity of the 
materials was admitted to be 4.1, 3.0 and 80.0, respectively 
for the acrylic pipe, oil and water contained within the pipe. 


During the process of image reconstruction, a structured 
mesh was used, subdividing the region under study in two 
subregions — inner portion and pipe wall as can be seen in 
Figure la) and 1b). In order to avoid the “inverse crime” 
described in Wirgin (2004), we used a more refined mesh to 
compute the simulated capacitance values and a less refined 
mesh for the solution of the inverse problem. 


Fig.l. Discretization domain with (a) 128 elements and (b) 
32 elements inside the pipe. 


For the execution of the numerical routines, it was used 
a computer with an Intel® Core™ i7-4790, 16 GB of RAM, 
running Ubuntu operational system and GNU/Linux 64 bits 
kernel 15.10. 


A cluster was assembled using a set of virtual machines 
on VirtualBox 4.1.8 running OpenMP protocol library. 


The processing algorithm was coded in Fortran 90 
language, using the GNU Fortran compiler. For pre- and 
post-processing, routines developed in Scilab language 
have been implemented. 
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Evaluation criteria 


Two quality assessment methods were employed to 
evaluate the obtained images: qualitative and quantitative. 
The qualitative assessment involves visual comparisons of 
similarity for different reconstruction methods. 


Reconstructed images obtained from the pipe section are 
also compared quantitatively using the following 
estimators: 


a) Normalized mean square error (NMSE): 


2 
N ref_ rec 
_ AG —ēi ) 


Ae = — 7 
mi (ear m 
b) Normalized absolute error — NAE: 
zh ere ere 
ôe = ——_— (8) 
re” | 
c) Correlation coefficient — Rw: 
yu rec_zrec)( ,tef_sref 
Re — i AG E JG E ) (9) 


1 
(Eer -eree P) (EN (e-e) 

Where eff and e7® are respectively the permittivity 
value of the element i for the reference dielectric 
distribution and reconstruted dielectric distribution by LBP 
algorithm and PSO, and é™ and é' are their average 
values, respectively. 


The best algorithm is the one with small values of ĝe 
and Ae, and values of Rx, close to unity. The normalized 
quadratic error is sensitive to large errors of some elements, 
while the normalized absolute error is sensitive to small 
errors on many factors, where the correlation coefficient 
indicates the spatial similarity between the reference image 
and the reconstructed image. 


IHI. RESULTS 


Images of the simulated flows profiles are shown on 
Figure 2. Those images should be compared with Figures 3, 
4 and 5 for the qualitative evaluation of the results. 
Numerical noise of 3% and 5% were also added to the 
capacitance values in order to test the robustness of 
reconstruction algorithms. 


In order to compare the quality of images obtained by 
the PSO algorithm, we employed linear back projection 
(LBP) algorithm as reference for the simulated flows shown 
on Figure 2. Based on the sensitivity matrix model, LBP is 
still the most cited reconstruction method and in its simplest 
implementation it presumes that the sensitivity matrix is 
invariant within the studied area (Li, 2015). All images were 
accomplished using a 5000 particles swarm and 500 
iterations. 
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Numerical evaluation for the chosen estimators is listed 
on Tables 1, 2 and 3 for capacitance values without 
numerical noise, 3% and 5% of noise, respectively. 


(a) (b) (c) 


Fig.2. Dielectric simulated flow profiles: (a) core flow, (b) 
bubbles and (c) stratified. 


Table 1. Error values e related to the image reconstruction 
process without numerical noise. 


Flow Metho | NMSE | NAE P 
Pattern d (%) (%) a 
LBP 57.90 | 56.97 | 92.96 
1 1 7 
CORE 
PSO 0.137 | 0.112 | 99.99 
9 
LBP 103.5 | 125.6 | 22.89 
10 70 0 
BUBBLE 
PSO 0.010 | 0.010 | 100.0 
00 
LBP 51.75 | 32.23 | 85.89 
STRATIF 8 0 1 
IED PSO 2.842 | 1.883 | 99.96 
1 


Table 2. Error values € related to the image reconstruction 
process with numerical noise of 3%. 


Flow Metho NMSE NAE 


Pattern d (%) (%) Ry 
LBP 59.134 | 58.560 91.836 
CORE 
PSO 9.831 4.472 | 99.548 
LBP 102.462 | 123.595 | 22.653 
BUBBLE 


PSO 36.812 | 24.552 | 93.523 
STRATIFI LBP 51.189 | 31.665 | 86.197 
ED PSO 6.281 3.400 | 99.815 
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Table 3. Error values £ related to the image reconstruction 
process with numerical noise of 5%. 


Flow Metho | NMSE NAE 


Ry 
Pattern d (%) (%) 
LBP 62.440 | 64.772 | 90.778 
CORE 
PSO 18.389 7.613 | 98.400 
LBP | 102.998 | 122.218 | 21.327 
BUBBLE 


PSO 55.840 37.547 | 85.370 


STRATIFI | LBP | 52.236 | 33.208 | 85.652 


ED PSO 9.856 5.542 | 99.547 


4 e) ~ 
(b) (c) 

» A i 
(e) (f) 


Figure 3. Reconstructed images from simulated values 
without numerical noise. 


ý & ~ 
(b) (c) 

tf eal i 
(e) (f) 


Figure 4. Reconstructed images from simulated values 
without numerical noise of 3%. 


O ý% A ~ 
(a) (b) (c) 
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ns ey al k 
(d) (e) (f) 


Fig.5. Reconstructed images from simulated values 
without numerical noise of 5%. 


IV. DISCUSSION 


PSO technique applied to ECT imaging has compelling 
results concerning its potentiality. Even in the presence of 
noise, it is possible to observe, not only visually but also by 
the quality indicators above, a good agreement between the 
proposed permittivity models and the accomplished results 
for the distinct flow patterns tested. 


vV. CONCLUSION 


However, one of the main limitations of PSO technique 
is reached when the complexity of the optimization 
functional increases in conjunction with the amount of 
image pixels and, therefore, the number of parameters to be 
determined. This situation limits the maximum image 
resolution that can be achieved in real time applications. A 
possible solution is to grow the swarm size, which facilitates 
the determination of an optimal solution, but such approach 
is time consuming when a single processing unit is 
employed. 
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